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Distributed cortical solutions of magnetoencephalography (MEG) and electroencephalog- 
raphy (EEG) exhibit complex spatial and temporal dynamics. The extraction of patterns of 
interest and dynamic features from these cortical signals has so far relied on the expertise 
of investigators. There is a definite need in both clinical and neuroscience research for a 
method that will extract critical features from high-dimensional neuroimaging data in an 
automatic fashion. We have previously demonstrated the use of optical flow techniques 
for evaluating the kinematic properties of motion field projected on non-flat manifolds like 
in a cortical surface. We have further extended this framework to automatically detect 
features in the optical flow vector field by using the modified and extended 2-Riemannian 
Helmholtz-Hodge decomposition (HHD). Here, we applied these mathematical models on 
simulation and MEG data recorded from a healthy individual during a somatosensory exper- 
iment and an epilepsy pediatric patient during sleep. We tested whether our technique can 
automatically extract salient dynamical features of cortical activity. Simulation results indi- 
cated that we can precisely reproduce the simulated cortical dynamics with HHD; encode 
them in sparse features and represent the propagation of brain activity between distinct 
cortical areas. Using HHD, we decoded the somatosensory N20 component into two HHD 
features and represented the dynamics of brain activity as a traveling source between 
two primary somatosensory regions. In the epilepsy patient, we displayed the propagation 
of the epileptic activity around the margins of a brain lesion. Our findings indicate that 
HHD measures computed from cortical dynamics can: (i) quantitatively access the corti- 
cal dynamics in both healthy and disease brain in terms of sparse features and dynamic 
brain activity propagation between distinct cortical areas, and (ii) facilitate a reproducible, 
automated analysis of experimental and clinical MEG/EEG source imaging data. 



Keywords: motion field, optical flow, MEG source imaging, Helmholtz-Hodge decomposition, epilepsy 



1. INTRODUCTION 

MEG and EEG are the most direct correlates of neural cur- 
rents measured externally (Baillet et al., 2001). Recent advances 
in both hardware and software (i.e., increase in number of 
sensors, faster microprocessors, and more accurate cortical sur- 
face reconstructions) have led to significant enhancement in the 
temporal and spatial resolution of these methods, which can 
now reach, sub-millisecond and sub-centimeter levels respec- 
tively (Murray et al, 2008; Papadelis et al., 2009). MEG and 
EEG measure the magnetic and electric correlates of intra- 
cranial currents respectively. In order to estimate the location 
and time-course of these neural current generators, we need 
to solve an ill-posed and non-unique inverse problem. The 
non-uniqueness of the inverse problem is a result of the non- 
triviality in the quasi-static Helmholtz equation that links the 



intra-cranial current sources to the observed extra-cranial fields. 
Spatiotemporal distributed source solutions, like minimum-norm 
estimate (MNE), of MEG/EEG have been proposed to over- 
come this non-triviality (Dale and Sereno, 1993; Hamalainen 
and Ilmoniemi, 1994). The resulting current distribution incor- 
porates the anatomical information for each individual brain 
from the magnetic resonance imaging (MRI), and calculates 
the time-course of source distributions usually constrained to 
the cortex. These MNE solvers leads to spatiotemporal linear 
solutions, and have been extensively used in the neuroimaging 
community for their relatively accurate source localizations and 
robustness to the noise levels normally present in MEG/EEG 
data sets. 

This type of analysis leads to a huge amount of high- 
dimensional data containing large information in both time and 
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space. Current approaches normally relies on studying cortical 
current variations at selected short latencies or by subtracting 
experimental conditions (e.g., standard minus deviant) to find fea- 
tures of interest in both space and time. These approaches allow 
mapping the local traveling of spatiotemporal cortical current 
activations on the cortical manifold. This propagation of brain 
activity via surface connections may represent propagating waves 
of cortical activity, which can emerge, transverse, or contract on 
the cortical surface (Ermentrout and Kleinfeld, 2001; Roxin et al., 
2005; Gramfort etal.,2011). 

The extraction of salient features mapping the spatiotemporal 
propagation of brain activity across different cortical regions relies 
so far on the expertise of neuroscience investigators or clinicians, 
who visually identify and quantify these patterns by using sta- 
tistical tools. However, this procedure remains problematic since 
the statistical models are prone to subjective bias of the investiga- 
tor. There is a definite need for methods that allow the automatic 
and a priori extraction of features of interest in evolving cortical 
dynamics and mapping in an automatic fashion the propagating 
activity between different cortical regions. 

In this paper, we propose a novel method to overcome these 
difficulties using mathematical formalism we described previously 
(Lefevre and Baillet, 2008; Khan et al, 2011). Our method extracts 
the spatiotemporal dynamics of EEG/MEG cortical sources using a 
combination of 2-Riemannian optical flow and Helmholtz-Hodge 
decomposition (HHD) on a highly-curved cortical manifold. The 
optical flow is a computer vision technique that represents the 
apparent motion in the time series of images. The HHD can 
automatically extract salient features from the optical flow. 

Mathematically; HDD decomposes optical flow into: 

• a non-rotational component deriving from the gradient of a 
scalar potential U; 

• a non-diverging component deriving from the rotational of a 
scalar potential A (resp. vectorial potential, in 3D); and 

• an harmonic part H, i.e., whose Laplacian vanishes. 

In HHD, formalism features of interest are represented as crit- 
ical points of scalar field U and A. Finding features as critical 
points in global field potential is much less sensitive to noise 
in the data and therefore less likely to get false positives (Khan, 
2010). In comparison to current density, its spatiotemporal diver- 
gent component U yields more focal and compact representation 
of the cortical activity (Slater et al, 2008; Khan et al, 2011). In 
order to estimate the spatiotemporal propagation of brain activ- 
ity, we should initially estimate and extract the features of interest 
using the diverging component U. Subsequently, HHD's harmonic 
part H can infer how the information propagates between cortical 
areas, by a vector field which is both irrotational and incompress- 
ible. This Laplacian vector field can explain causal effects exerted 
by one brain region onto another. Particularly this vector field is 
especially applicable in revealing dynamics, which occur briskly in 
time and over short distances on the cortical manifold. 

The detection of features in optical flow motion field using 
HHD has already been applied in many different imaging fields 
(Palit et al, 2005; Guo et al, 2006). For instance, it is used in car- 
diac video analysis, to detect features in cardiac motion fields that 



reflect pathological activity in the dynamics of cardiac electrical 
activity. The proposed method is based on previous work done by 
our group (Lefevre and Baillet, 2008) where we introduced a vari- 
ational method to estimate the optical flow on non-flat surfaces 
using a Riemannian formulation. 

This previously proposed technique was used to analyze the 
global dynamics of cortically-distributed source images obtained 
from MEG or EEG data with also limited quantification of local 
dynamics (Lefevre and Baillet, 2009). It was recently extended by 
introducing a new formalism to detect local features of the optical 
flow of cortical dynamics using a modified and extended approach 
to HHD (Slater et al, 2008; Khan et aL, 2011). 

This paper is structured as following: the Riemannian frame- 
work for optical flow on non-flat surfaces is first briefly introduced. 
The HHD formalism on 2-Riemannian manifolds will be pre- 
sented next. Lastly, the application of HHD on simulated and 
human MEG data will be presented. The methods discussed in 
this paper are implemented in MatLab and are available for down- 
load as plugin to Brainstorm (MEG/EEG data processing software) 
(Tadel et al., 2011) at http://neuroimage.usc.edu/brainstorm. 
These methods wUl also soon be available for MNE-Python 
framework (Gramfort et al., 2014). 

2. MATERIALS AND METHODS 

The HHD-based feature detection technique consists of three 
distinct steps. 

• First the optical flow of distributed MEG/EEG MNE estimates 
is computed on the cortical manifold. In Section 2.1, we will 
briefly introduce optical flow and its mathematical formulation. 

• In the second step, HHD is applied on optical flow computed 
previously. We will present HHD framework in Section 2.2 and 
concisely describe its axioms. 

• Lastly, detecting the feature of interest in cortical dynamics now 
becomes the simple problem of identifying critical points in 
HHD scalar potential U. The traveling cortical dynamics can be 
tracked by vectors having highest norm in vector field H (see 
Section 2.3 for details). 

2.1. OPTICAL FLOW 

We have introduced the concept of optical flow on a 2-Riemannian 
manifold (Lefevre and Baillet, 2008), and we shall briefly summa- 
rize the approach as follows. Under the seminal hypothesis of the 
conservation of a scalar field 1 along streamlines, defined on a 
surface Ad, the optical flow V is a vector field that satisfies: 

d,I + g{\,VMl) = Q- (1) 

Note that the scalar product g(.,.) is modified by the local curva- 
ture of AA-, the domain of interest. The solution to Equation (1) 
is not unique as long as the components of V(p, f ) orthogonal 
to V M I arc left unconstrained. This so-called "aperture problem" 
has been addressed by a large number of methods using e.g., regu- 
larization approaches. These latter approaches may be formalized 
as the minimization problem of an energy functional, which both 
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includes the regularity of the solution and the agreement to the 
model: 



f(V) 



/ 

Jm 



di 
9f 



+ ^(V,Va4/) 



dii+X f C(V)d/x. (2) 

Jm 



Here, we considered a regularity factor operating quadratically on 
the gradient of the expected vector field: 



C(V) = rr('VV.VV). 



(3) 



Note that in order to be an intrinsic tensor, the gradient of a vec- 
tor field must be defined as the covariant derivative associated to 
the manifold A4. Due to space constraints, we need to refer to 
essential elements of differential geometry for more information 
on this notion (Do Carmo, 1993). 

2.2. HELMHOLTZ-HODGE DECOMPOSITION ON 2-RIEMANNIAN 
MANIFOLD 

Let us consider 7V4 as a surface (or manifold) parameterized by 
local charts (xi, X2). It is thus possible to obtain a normal vector 
at each surface point: 



dxi 8x2 

Note that the normal vector does not depend on the choice of 
the parameterization {xi, xj). We then define the gradient and 
divergence operators through duality: 

d[/(V) = g(VA^U,V), 

\ UdwMii = -[ g(n,WMU). 
JM JM 

Finally, scalar and vectorial curl operators are obtained through: 

CumA = V^kA a n, 
cu^ H = div^K (H A n) . 

Again, these formulas are intrinsic expressions that do not depend 
on the parameterization of the surface. 

Let us then reformulate results established in Polthier and 
Preuss (2003) for Riemannian manifolds. Given V, a vector field 
in r^(M), there exists unique functions U and A in L^(A4) and 
a vector field H in T ^ {M ) such that: 



V = Va^[/ + Cua^A + H, (4) 

CUA^ (Va^ L^) = 0, divA^ (CuA^ A) = 0, 
divA^H = 0, CUA^H = 0. 



with 



Following classical constructions, U and A minimize the two 
following functional: 

[ \\y-^MU\\\ [ ||V-CUA^A||2, 
JM JM 



where 11-11 is the norm associated to the Riemannian metric gi.,.). 
These two functionals are convex and therefore have a minimum 
on L^{A4) satisfying: 

V,peLHM), ( giY,VM<P)=i g(^MU,^M<P), (5) 

JM JM 



y<l>eL'{M), [ ,g(V,CuAi0) 
JM 



JM 



^(CuAlA,CuAi</)). 



(6) 



Through (0i, .. .,</>„) as basis functions, we may write 
\J=iUi,. . .,U„)'^,A=iAi,. . .,A„f, and Equations (5) and (6) 
read, using array notations: 



/ g('^M(t>i,^M4>j) 
JM 

/ ^(CuA40,,CuA40j) 
JM 



U : 



[ g(V,VA40O 
JM 

[ g(V,CviM<P,) 
JM 



(7) 



(8) 



The harmonic component H of vector field V is simply obtained 
through: 



H = V - Vai U - CuA^A. 



(9) 



2.3. FEATURE DETECTION AS CRITICAL POINTS OF HHD POTENTIALS 

The critical points of a vector field are often classified depending 
on the eigenvalues of the Jacobian matrix defined locally in a vec- 
tor field. In our case, however, critical points of the flow can be 
found as local extrema of the potentials A and U. Finding features 
as critical points on global potential fields is much less sensitive 
to noise in the data and therefore is less susceptible to false pos- 
itives, than with methods using local Jacobian eigenvalues (Tong 
et al., 2003). Moreover unlike eigenvalues methods, HHD do not 
pre-specify the number of features. 

A sink (resp. a source) corresponds to a local minimum (resp. 
maximum) of U. Similarly, a counter-clockwise (resp. clockwise) 
vortex may be detected through a local minimum (resp. maxi- 
mum) of A. Detection of traveling cortical activity on the surface 
can be performed by tracking vector H having largest norm. 

3. RESULTS 

We wiU now present the application of HHD on one simulated 
and two actual MEG datasets. In Section 3.1 under a simulated 
scenario, HHD is used to track and encode cortical activity as it 
emerges from the somatosensory cortex, traveling along the cen- 
tral sulcus, and receding in the inferior frontal gyrus. The first 
MEG dataset presented in Section 3.2 is obtained from a median 
nerve stimulation paradigm that consists of a train of electrical 
stimuli applied on the wrist of a healthy adult individual. It is a 
typical experimental paradigm to elicit activity within the primary 
somatosensory cortex. The second MEG dataset shown in Section 
3 .3 is from a pediatric epilepsy patient with tuberous sclerosis com- 
plex. In this dataset, we track the propagation of the epileptogenic 
activity within the irritative zone. 

3.1. SIMULATION DATA 

The use of HHD to address the quantitative and qualitative 
evaluation of this technique is illustrated below. 
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3. 1. 1. Generation 

The brain surface from Freesurfer's FsAverage was selected to 
demonstrate the apphcabihty of the method on the cortical mani- 
fold. This surface consisted of 10,242 vertices and 20,480 triangles. 
A source on this manifold was generated in the vicinity of the cen- 
tral sulcus, also known as the primary somatosensory area (SI). 
This source was grown to a patch of 5 cm^ (geodesic area) in 
five time steps. Subsequently, a constant vector field was defined 
in time from the vectorial heat equation. An advection equation 
(Lefevre and Baillet, 2008) was used to transverse this patch on the 
manifold. Finally, this patch was contracted in five time steps in the 
vicinity of the inferior frontal gyrus. The results of this simulation 
are presented in Figure 1. The three stages of the source evolution 
are shown as Figures lA-C. 

3.1.2. Analysis 

We first applied the optical flow on this simulated activity. Opti- 
cal flow transformed the dynamics of the source's evolution in 
terms of the motion vector fields that were emerging, traversing, 
and receding. This optical flow computation was performed in the 




FIGURE 1 I Examples of different types of vector fields (green arrows) 
and their corresponding U and (H) HIHD components. (A) Emerging 
source represented by a minima in the scalar potential U. (B) Traveling 
source detected by (H). (C) Receding source represented by a maxima in 
the scalar potential U. 



MNE-Python framework and took 10.75 s to compute. We then 
applied HHD to the optical flow at each time step, which detects 
three main features corresponding to the three stages of the simu- 
lation. HHD took 30.36 s to decompose optical flow in the discrete 
feature set. The simulation was performed on a workstation having 
a Dual Octa Core Xeon CPUs (32 Threads) with 64 GB of RAM. 

3.1.3. Results 

Figure 1 shows the applicability of HHD on the distributed cur- 
rent density maps for three different source configurations: (A) an 
emerging source; (B) a moving source having constant velocity; 
and (C) a descending source (sink). HHD automatically detects 
these three features from optical flow as critical points in U and H. 
A source represented by a minima in the scalar potential U (texture 
colormap in blue) is shown in Figure lA. The vector field (green 
arrows) represents optical flow computed earlier. It is indicated in 
Figure lA that HHD is able to capture the growing dynamics of 
the cortical activity. In Figure IB, traveling activity is tracked by 
the highest norm vectors in field H. The texture colormap repre- 
sents the norm of H, where the highest norm is shown with the 
dark violet texture map. In Figure IC, receding cortical dynamics 
are captured by a maxima in the scalar potential U. This texture 
colormap represents U with a maxima shown in red. Again, the 
vector field is shown with green arrows representing optical flow. 

3.2. FEATURE ANALYSIS OF EXPERIMENTAL MEG SOMATOSENSORY 

DATA 
3.2.1. Experiment 

MEG data was recorded from a 28-year-old healthy female individ- 
ual at the MEG laboratory of the Center for Mind/Brain Sciences 
(CIMeC), University of Trento, Italy. During the experiment, the 
median nerves of her right and left wrists were electrically stimu- 
lated transcutaneously. Approval was obtained from the University 
of Trento Ethics Committee, Italy, and the subject gave her writ- 
ten informed consent before the experiment. Constant current 
stimuli with a duration of 0.2 ms, and a pseudo-randomized 
inter- stimulus interval of 250 ± 50 ms were applied to the sub- 
ject. Before the start of the experiment, we measured two basic 
intensity thresholds for each of the subject's wrists. The two para- 
meters were the motor threshold (MTH), defined as the minimal 
stimulus intensity needed to produce thumb twitching, and the 
sensory threshold (STH), defined as the minimal stimulus inten- 
sity at which the subject was just able to feel a train of stimulus 
pulses. The stimuli were delivered either to the right or to the 
left wrist with two intensity levels of M = MTH + 0.25 A and 
S = STH + 0.25 A, where A = MTH - STH. Here, we used only 
the data from the right median nerve stimulation and the M inten- 
sity level. More details about the experimental design can be found 
inPapadelis et al. (2012). 

Somatosensory evoked fields (SEFs) were recorded at a sam- 
pling rate of 5 kHz by using a 306-channel (204 first order pla- 
nar gradiometers, 102 magnetometers) VectorView MEG system 
(Elekta-Neuromag Ltd., Helsinki, Finland) placed inside a magnet- 
ically shielded room (AK3B, Vacuumschmelze, Hanau, Germany). 
Hardware filters were adjusted to band-pass the MEG signal in the 
frequency range of 0.01-1000 Hz. Data from 200 trials were used 
in this study. 
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3.2.2. Data analysis 

To compensate for head movements during the measurements and 
suppressing external magnetic disturbances, the signal space sepa- 
ration (SSS) algorithm (Taulu and Simola, 2006), as implemented 
with the MaxFilter software (Elekta-Neuromag), was performed 
offline on the raw MEG data. The corrected MEG data were 
then filtered offline in the band of 0.1-150 Hz and epoched from 
—50 to 200 ms relative to the stimulus onset. Trials contaminated 
with artifacts were excluded from further processing. In total, 160 
trials survived the rejection criteria. These trials were baseline cor- 
rected and then averaged. The subject's MRI was processed using 
FreeSurfer; cortical and head surfaces were extracted. Cortical sur- 
face was downgraded to 50,003 points and elementary current 
dipoles were positioned at the surface of the cortex of the subject. 
The multi-sphere forward model was computed and the standard 
minimum-norm was used to estimate cortical currents. The pre- 
processing and generation of both the forward and inverse solution 
was done in Brainstorm (Tadel et al., 2011). The MNE solution 
was computed at every point within the time window that repre- 
sents the N20 component (from 15 to 25 ms) (see Figure 2, upper 
panel). The N20 component represents the first cortical response 



to the electrical stimulation of the median nerve. Optical flow 
was estimated from the minimum-norm data. HHD was then 
applied on this optical flow to detect sparse features in cortical 
dynamics. The computation of optical flow took 50.33 s whereas 
HHD took 123.73 s on the workstation mentioned in Section 

3.1.2. This application was done in the MatLab implementation 
of the HHD. 

3.2.3. Results 

Figure 2 shows the applicability of HHD on the distributed source 
maps of MEG data recorded during the median nerve stimulation 
experiment. Our method reveals the diverging and contracting 
cortical mechanisms in the primary somatosensory cortex. Brain 
responses during N20 were automatically decomposed into two 
features: a source and a sink. 

In Figure 2, the upper panel shows the butterfly plot of 
somatosensory evoked fields (SEFs). Two features identified by 
HHD are shown in the middle panel of Figure 2; an emerging 
source at 17.3 ms (MNI: -46.74, -30.17, 66.84) and a descending 
source (sink) at 23.2 ms (MNI: -55.86, -22.26, 51.20). The green 
arrows represent the motion field of the cortical dynamics as 




FIGURE 2 I Decomposition of cortical activity in a series of dynamic features. Response to median nerve stimulation of tine right wrist. Upper panel shows 
the butterfly plot of somatosensory evoked fields (SEFs). Two features identified by HHD are shown in the middle panel. Bottom panel shows zoom view of the 
detected features. 
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estimated by the optical flow. HHD reveals that the first corti- 
cal responses appear at the top of SI in the central sulcus, travels 
down the sulcus, and then sinks down at lower edge of the sul- 
cus. The bottom panel presents a zoom up view of these two 
features. HHD reduces the size of the dataset from 50,003 x 50 
(spatial dimension x temporal dimension) to two main features. 
This example presents the concept of our methodology to obtain 
a compact representation of cortical activity during a cognitive 
experiment. 

3.3. HHD CHARACTERIZATION OF EPILEPTIC ACTIVITY 

Our second application of HHD is on MEG data acquired from a 
subject having clinical history of refractory epilepsy. In this appli- 
cation, we will focus on the property of the [/ component of HHD, 
which allows for the detection of activity that changes rapidly in 
both time and space. 

3.3. 1. Epileptic data 

The data was recorded from a 20-month-old boy, who presented 
his first seizure at the age of 3 months, with refractory epilepsy 
as a result of TSC. The patient had multiple subcortical tubers 
identified on his MRI as patchy areas of T2 prolongation, sta- 
ble over time (see Figure 3). Long-term monitoring revealed 54 
seizures in total (46 electroclinical and 8 electrographic seizures) 
with duration of 10-44 s. The seizure onset was localized at the 
right posterior quadrant (electrodes P8, 02, P4, T8, and C4). Rou- 
tine and ambulatory EEG has indicated frequent interictal sharp 
waves at electrodes C4, Pz, P4, and P8, as well as slowing at the 
right posterior quadrant. 

MEG data were recorded for 45 min during sleep at the 
BabyMEG facility located at the Radiology Suite of Boston's 
Children Hospital (Waltham, MA, USA). MEG recordings were 
performed using a 74-sensor MEG system especially designed 
for pediatric use ("babySQUID" - Tristan Technologies Inc., 
San Diego, CA, USA). The babySQUID system is accommo- 
dated in a single-layer magnetically shielded room (MSR). MEG 
data was sampled at 1024 samples per second. The sensor 
array of MEG was covering the right posterior quadrant. Tl- 
weighted high-resolution magnetization-prepared rapid gradient 
echo (MPRAGE) structural images were acquired on a 3.0-T 
Siemens Trio whole body MR scanner (Siemens Medical Systems, 
Erlangen, Germany) using a 32 channel head coil. Details about the 
experimental procedure can be found elsewhere (PapadeHs et al., 




FIGURE 3 I Lesion is sliown onT2 (left) andT1 (righit) withi red outline. 

Epileptic activity emerging from the edge of the lesion is shown as texture 
map. 



2013). Research MEG and MRI data were acquired and analyzed 
after explicit parental consent under a protocol approval by the 
Boston's Children Hospital institutional review board. 

3.3.2. Analysis 

A high number of interictal spikes (>20) were identified by a 
board-certified epileptologist with a consistent spatiotemporal 
pattern indicating a focal source in the right posterior quadrant. 
The MEG data was then filtered offline in the band of 0. 1-145 Hz. 
The subject's MRI was processed using brainvisa; cortical and 
head surfaces were extracted. To compute the forward solution, 
a boundary-element model (BEM) with a single compartment 
bounded by the skull's inner surface was assumed (Hamalainen 
et al, 1993). The watershed algorithm was used to generate the 
inner skull surface triangulations from the high-resolution Tl MR 
images of each participant. The current distribution was estimated 
using the MNE by fixing the source orientation to be perpendicu- 
lar to the cortex. The noise covariance matrix used to calculate the 
inverse operator was estimated from empty-room data. In order 
to reduce the bias of the MNEs toward superficial currents, we 
incorporated depth weighting by adjusting the source covariance 
matrix (Lin et al., 2006). To estimate the epileptic foci for the 
subject from HHD, a single sharp epileptic spike (Figure 4A) was 
selected. MNE was then computed for this spike in both volume 
and cortical space. Optical flow was calculated from the MNE esti- 
mated cortical currents at each time point during a time window 
shown in Figure 4A. HHD was computed by decomposing this 
optical flow and extracting the diverging scalar field L^. The com- 
putation of optical flow and HHD took 34.7 s on the workstation 
mentioned in Section 3.1. 

3.3.3. Results 

In Figure 4B, MNE activity is presented in the volume at the two 
arrow points in Figure 4A. Average MNE activity on the corti- 
cal manifold at these two points is shown in Figure 4C. Critical 
points are then searched in the scalar field. This process results 
in finding two critical points in this time window, which corre- 
sponds to the signature features of this epileptic spike. A source in 
the posterior occipital region SI, from where the activity emerges, 
is shown in blue in Figure 4D. This activity sinks in the anterior 
parietal cortex in the vicinity of a subcortical tuber. This sinking 
activity is shown in red in Figure 4D. Finally, the diverging vec- 
tor field Vjvf U is computed by taking the gradient of the scalar 
field U. In Figure 4D, the divergence vector field for this epileptic 
spike is shown with black arrows whereas color texture represents 
strength of U. 

4. DISCUSSION 

In this paper, we present the applicability of HHD to high- 
dimensional neuroimaging data. The HHD-based sparse feature 
encoding technique works in three steps. First, MNE is computed 
to estimate cortical current activity. Optical flow is then used to 
estimate the motion field of distributed cortical dynamics. Finally, 
the optical flow is decomposed into sparse and compact features 
using HHD and the neuroimaging feature extraction simplifies to 
the problem of finding critical points in the scalar potential and 
by the highest norm vector of H. 
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FIGURE 4 I Encoding of epileptic spike in diverging features. (A) W 

traces witli epileptic spike marked. (B) MRI with MEG activity 
represented as the color texture. (C) Average MEG activity during the 
spike. (D) U HHD scalar potential is mapped onto the cortical surface 



nAm 



EG using textured colors. The divergence part of vector flow of MEG sources 

is represented by green arrows at each vertex location. Critical points in 
the U map (shown with magenta sphere) reveal sources shown in dark 
blue and sink in red. 



The method has been optimized for running in muhi-core 
CPUs in order to decompose high-dimensional MEG/EEG data 
(~GB) in few seconds on a highly dense (-50,000 vertices) corti- 
cal mesh. The method provides a compact representation of the 
cortical dynamics. The method also allows to track the cortical 
activity as it appears, disappears, or travels on the cortical mani- 
fold. Moreover, HHD can also infer the information flow between 
cortical regions over short distances and times. Here, we present 
the appHcation of our method on both simulated and actual MEG 
data sets. 

In the simulated scenario, we generated synthetic cortical activ- 
ity from the heat vector field and the advection equation. We then 
applied optical flow and HHD on this simulated data compacting 
the different stages of these dynamics in three sparse features. 
These results were in absolute accordance with our generated 
data. In the cognitive neuroscience appKcation, we encoded spa- 
tiotemporal time varying maps of cortical activity during electrical 
stimulation of the median nerve. We showed that HHD and opti- 
cal flow can compress cortical dynamics of N20 component in 
two distinct features; a source and a sink. This is in accordance 
to recent evidence in the literature of somatosensory processing 
which demonstrate that the earliest cortical response component 



(N20) after median nerve stimulation is generated by two gener- 
ators located within SI subdivisions and being active following 
a serial fashion (Inui et al., 2004; Papadelis et al, 2011). The 
application of our method in neuroscience data may enable the 
reproducibility of cognitive experiments across different sessions 
or research centers. 

Our method can significantly contribute to epilepsy research, 
because it is able to detect and map the propagation of interictal 
spikes across time. MEG and EEG studies have so far shown that 
these two neuroimaging methods can non-invasively detect the 
propagation of spikes in epilepsy patients (Sutherling and Barth, 
1989; Emerson et al, 1995; Tanaka et al., 2010). We also like to 
emphasize the importance of initial detection of the spike with a 
consistent spatiotemporal pattern indicating a focal source by an 
epileptologist, as the method seeds from it. 

The investigation of epileptic spike propagation is important 
for the understanding of the pathophysiology of epilepsy and for 
the appropriate clinical decision during the presurgical evaluation 
of epileptic patients. Spike propagation can reflect neural networks 
associated with epilepsy (Spencer, 2002), while the propagation 
pattern of interictal spikes has been shown to be related to the 
outcome of epilepsy surgery (Hufnagel et al., 2000; Schulz et al. 
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2000). In the application of HHD to epilepsy data, we exploited 
the sensitivity of the divergence component U to fast changes 
in space and time, such as these during a propagating interictal 
spike. The epileptic activity at the onset of the interictal spike 
can be sharp spatio-temporarily, and may easily be detected by 
HHD. We used HHD in deciphering these complex interictal spike 
dynamics. Using HHD, we identified where the spike initiated but 
also tracked its activity as it propagated on the cortical manifold. 
The application of HHD found two critical points on the corti- 
cal manifold, a source and a sink, both located at the abnormally 
developed tissue surrounding a tuber rather than the tuber itself 
The epileptogenic activity was propagating across time along the 
borders of the tuber, and in any instance was crossing the tuber 
itself. Our findings are in agreement with previous studies indicat- 
ing that in TSC patients with epilepsy the epileptogenic tissue is 
predominately localized in the surroundings of the cortical tubers 
(Weiner, 2004; Xiao et al., 2006; Major et al, 2009), and a single 
case study published in this issue (Hunold et al., 2014). The high 
sensitivity of our method allowed us to map the evolution of the 
epileptiform activity across time with respect to the location of the 
tuber. This critical, automatically extracted, spatiotemporal infor- 
mation of interictal spikes may provide more accurate information 
of spike propagation in epilepsy patients compared to the classical, 
observer-dependent methods, and thus it may be clinically useful 
in the presurgical evaluation of epilepsy patients. 

5. CONCLUSION 

We demonstrate the applicability of our HHD technique on high- 
dimensional electrophysiological data from neuroscience and clin- 
ical research data. Salient features of our technique, which are 
demonstrated by our results are: (i) sensitivity to spatiotempo- 
ral diverging cortical sources rapidly evolving in space and time 
(within few milliseconds), (ii) consideration of geometry of the 
cortical manifold on which the neural activity is evolving, (iii) 
automatic extraction of the spatiotemporal features, (iv) auto- 
matic characterization of the cortical activity propagation across 
different brain regions, (vi) visualization of the salient feature of 
cortical activity, and (vii) application in both cognitive and clini- 
cal neuroscience (i.e., propagation of epileptiform activity). Here, 
we present some of the possible applications of HHD in the neu- 
roimaging field. We strongly believe that there are much more 
applications of our method in both neuroscience as well as clin- 
ical research. Further improvements of our method include: the 
discretization to higher-order finite element analysis, statistical 
analysis of the detected features, and the further reduction of the 
algorithm computational complexity. HHD is freely available as 
a plugin for major MEG processing suites (i.e.. Brainstorm and 
MNE-Python) and readers are encouraged to use and extend the 
proposed method. 
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